home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 34 / Amiga Format CD34 (1998-11-20)(Future Publishing)(GB)[!][Christmas issue].iso / -seriously_amiga- / programming / c / mesa-2.6 / src-glu / project.c < prev    next >
C/C++ Source or Header  |  1998-10-01  |  12KB  |  440 lines
  1. /* $Id: project.c,v 1.5 1997/07/24 01:28:44 brianp Exp $ */
  2.  
  3. /*
  4.  * Mesa 3-D graphics library
  5.  * Version:  2.4
  6.  * Copyright (C) 1995-1997  Brian Paul
  7.  *
  8.  * This library is free software; you can redistribute it and/or
  9.  * modify it under the terms of the GNU Library General Public
  10.  * License as published by the Free Software Foundation; either
  11.  * version 2 of the License, or (at your option) any later version.
  12.  *
  13.  * This library is distributed in the hope that it will be useful,
  14.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16.  * Library General Public License for more details.
  17.  *
  18.  * You should have received a copy of the GNU Library General Public
  19.  * License along with this library; if not, write to the Free
  20.  * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  */
  22.  
  23.  
  24. /*
  25.  * project.c
  26.  *
  27.  * Version 1.0  27 Jun 1998
  28.  * by Jarno van der Linden
  29.  * jarno@kcbbs.gen.nz
  30.  *
  31.  * File created from project.c ver 1.5 and glu.h ver 1.9 using GenProtos
  32.  *
  33.  */
  34.  
  35.  
  36. #ifdef PC_HEADER
  37. #include "all.h"
  38. #else
  39. #include <stdio.h>
  40. #include <string.h>
  41. #include <math.h>
  42. #include "gluP.h"
  43. #endif
  44.  
  45.  
  46. /*
  47.  * This code was contributed by Marc Buffat (buffat@mecaflu.ec-lyon.fr).
  48.  * Thanks Marc!!!
  49.  */
  50.  
  51.  
  52.  
  53. /* implementation de gluProject et gluUnproject */
  54. /* M. Buffat 17/2/95 */
  55.  
  56.  
  57.  
  58. /*
  59.  * Transform a point (column vector) by a 4x4 matrix.  I.e.  out = m * in
  60.  * Input:  m - the 4x4 matrix
  61.  *         in - the 4x1 vector
  62.  * Output:  out - the resulting 4x1 vector.
  63.  */
  64. static void transform_point( GLdouble out[4], const GLdouble m[16],
  65.                  const GLdouble in[4] )
  66. {
  67. #define M(row,col)  m[col*4+row]
  68.    out[0] = M(0,0) * in[0] + M(0,1) * in[1] + M(0,2) * in[2] + M(0,3) * in[3];
  69.    out[1] = M(1,0) * in[0] + M(1,1) * in[1] + M(1,2) * in[2] + M(1,3) * in[3];
  70.    out[2] = M(2,0) * in[0] + M(2,1) * in[1] + M(2,2) * in[2] + M(2,3) * in[3];
  71.    out[3] = M(3,0) * in[0] + M(3,1) * in[1] + M(3,2) * in[2] + M(3,3) * in[3];
  72. #undef M
  73. }
  74.  
  75.  
  76.  
  77.  
  78. /*
  79.  * Perform a 4x4 matrix multiplication  (product = a x b).
  80.  * Input:  a, b - matrices to multiply
  81.  * Output:  product - product of a and b
  82.  */
  83. static void matmul( GLdouble *product, const GLdouble *a, const GLdouble *b )
  84. {
  85.    /* This matmul was contributed by Thomas Malik */
  86.    GLdouble temp[16];
  87.    GLint i;
  88.  
  89. #define A(row,col)  a[(col<<2)+row]
  90. #define B(row,col)  b[(col<<2)+row]
  91. #define T(row,col)  temp[(col<<2)+row]
  92.  
  93.    /* i-te Zeile */
  94.    for (i = 0; i < 4; i++)
  95.      {
  96.     T(i, 0) = A(i, 0) * B(0, 0) + A(i, 1) * B(1, 0) + A(i, 2) * B(2, 0) + A(i, 3) * B(3, 0);
  97.     T(i, 1) = A(i, 0) * B(0, 1) + A(i, 1) * B(1, 1) + A(i, 2) * B(2, 1) + A(i, 3) * B(3, 1);
  98.     T(i, 2) = A(i, 0) * B(0, 2) + A(i, 1) * B(1, 2) + A(i, 2) * B(2, 2) + A(i, 3) * B(3, 2);
  99.     T(i, 3) = A(i, 0) * B(0, 3) + A(i, 1) * B(1, 3) + A(i, 2) * B(2, 3) + A(i, 3) * B(3, 3);
  100.      }
  101.  
  102. #undef A
  103. #undef B
  104. #undef T
  105.    MEMCPY( product, temp, 16*sizeof(GLdouble) );
  106. }
  107.  
  108.  
  109. static GLdouble Identity[16] = {
  110.    1.0, 0.0, 0.0, 0.0,
  111.    0.0, 1.0, 0.0, 0.0,
  112.    0.0, 0.0, 1.0, 0.0,
  113.    0.0, 0.0, 0.0, 1.0
  114. };
  115.  
  116.  
  117.  
  118. /*
  119.  * Compute the inverse of a 4x4 matrix.  Contributed by scotter@lafn.org
  120.  */
  121. static void invert_matrix_general( const GLdouble *m, GLdouble *out )
  122. {
  123.  
  124. /* NB. OpenGL Matrices are COLUMN major. */
  125. #define MAT(m,r,c) (m)[(c)*4+(r)]
  126.  
  127. /* Here's some shorthand converting standard (row,column) to index. */
  128. #define m11 MAT(m,0,0)
  129. #define m12 MAT(m,0,1)
  130. #define m13 MAT(m,0,2)
  131. #define m14 MAT(m,0,3)
  132. #define m21 MAT(m,1,0)
  133. #define m22 MAT(m,1,1)
  134. #define m23 MAT(m,1,2)
  135. #define m24 MAT(m,1,3)
  136. #define m31 MAT(m,2,0)
  137. #define m32 MAT(m,2,1)
  138. #define m33 MAT(m,2,2)
  139. #define m34 MAT(m,2,3)
  140. #define m41 MAT(m,3,0)
  141. #define m42 MAT(m,3,1)
  142. #define m43 MAT(m,3,2)
  143. #define m44 MAT(m,3,3)
  144.  
  145.    GLdouble det;
  146.    GLdouble d12, d13, d23, d24, d34, d41;
  147.    GLdouble tmp[16]; /* Allow out == in. */
  148.  
  149.    /* Inverse = adjoint / det. (See linear algebra texts.)*/
  150.  
  151.    /* pre-compute 2x2 dets for last two rows when computing */
  152.    /* cofactors of first two rows. */
  153.    d12 = (m31*m42-m41*m32);
  154.    d13 = (m31*m43-m41*m33);
  155.    d23 = (m32*m43-m42*m33);
  156.    d24 = (m32*m44-m42*m34);
  157.    d34 = (m33*m44-m43*m34);
  158.    d41 = (m34*m41-m44*m31);
  159.  
  160.    tmp[0] =  (m22 * d34 - m23 * d24 + m24 * d23);
  161.    tmp[1] = -(m21 * d34 + m23 * d41 + m24 * d13);
  162.    tmp[2] =  (m21 * d24 + m22 * d41 + m24 * d12);
  163.    tmp[3] = -(m21 * d23 - m22 * d13 + m23 * d12);
  164.  
  165.    /* Compute determinant as early as possible using these cofactors. */
  166.    det = m11 * tmp[0] + m12 * tmp[1] + m13 * tmp[2] + m14 * tmp[3];
  167.  
  168.    /* Run singularity test. */
  169.    if (det == 0.0) {
  170.       /* printf("invert_matrix: Warning: Singular matrix.\n"); */
  171.       MEMCPY( out, Identity, 16*sizeof(GLdouble) );
  172.    }
  173.    else {
  174.       GLdouble invDet = 1.0 / det;
  175.       /* Compute rest of inverse. */
  176.       tmp[0] *= invDet;
  177.       tmp[1] *= invDet;
  178.       tmp[2] *= invDet;
  179.       tmp[3] *= invDet;
  180.  
  181.       tmp[4] = -(m12 * d34 - m13 * d24 + m14 * d23) * invDet;
  182.       tmp[5] =  (m11 * d34 + m13 * d41 + m14 * d13) * invDet;
  183.       tmp[6] = -(m11 * d24 + m12 * d41 + m14 * d12) * invDet;
  184.       tmp[7] =  (m11 * d23 - m12 * d13 + m13 * d12) * invDet;
  185.  
  186.       /* Pre-compute 2x2 dets for first two rows when computing */
  187.       /* cofactors of last two rows. */
  188.       d12 = m11*m22-m21*m12;
  189.       d13 = m11*m23-m21*m13;
  190.       d23 = m12*m23-m22*m13;
  191.       d24 = m12*m24-m22*m14;
  192.       d34 = m13*m24-m23*m14;
  193.       d41 = m14*m21-m24*m11;
  194.  
  195.       tmp[8] =  (m42 * d34 - m43 * d24 + m44 * d23) * invDet;
  196.       tmp[9] = -(m41 * d34 + m43 * d41 + m44 * d13) * invDet;
  197.       tmp[10] =  (m41 * d24 + m42 * d41 + m44 * d12) * invDet;
  198.       tmp[11] = -(m41 * d23 - m42 * d13 + m43 * d12) * invDet;
  199.       tmp[12] = -(m32 * d34 - m33 * d24 + m34 * d23) * invDet;
  200.       tmp[13] =  (m31 * d34 + m33 * d41 + m34 * d13) * invDet;
  201.       tmp[14] = -(m31 * d24 + m32 * d41 + m34 * d12) * invDet;
  202.       tmp[15] =  (m31 * d23 - m32 * d13 + m33 * d12) * invDet;
  203.  
  204.       MEMCPY(out, tmp, 16*sizeof(GLdouble));
  205.    }
  206.  
  207. #undef m11
  208. #undef m12
  209. #undef m13
  210. #undef m14
  211. #undef m21
  212. #undef m22
  213. #undef m23
  214. #undef m24
  215. #undef m31
  216. #undef m32
  217. #undef m33
  218. #undef m34
  219. #undef m41
  220. #undef m42
  221. #undef m43
  222. #undef m44
  223. #undef MAT
  224. }
  225.  
  226.  
  227. /*
  228.  * Invert matrix m.  This algorithm contributed by Stephane Rehel
  229.  * <rehel@worldnet.fr>
  230.  */
  231. static void invert_matrix( const GLdouble *m, GLdouble *out )
  232. {
  233. /* NB. OpenGL Matrices are COLUMN major. */
  234. #define MAT(m,r,c) (m)[(c)*4+(r)]
  235.  
  236. /* Here's some shorthand converting standard (row,column) to index. */
  237. #define m11 MAT(m,0,0)
  238. #define m12 MAT(m,0,1)
  239. #define m13 MAT(m,0,2)
  240. #define m14 MAT(m,0,3)
  241. #define m21 MAT(m,1,0)
  242. #define m22 MAT(m,1,1)
  243. #define m23 MAT(m,1,2)
  244. #define m24 MAT(m,1,3)
  245. #define m31 MAT(m,2,0)
  246. #define m32 MAT(m,2,1)
  247. #define m33 MAT(m,2,2)
  248. #define m34 MAT(m,2,3)
  249. #define m41 MAT(m,3,0)
  250. #define m42 MAT(m,3,1)
  251. #define m43 MAT(m,3,2)
  252. #define m44 MAT(m,3,3)
  253.  
  254.    register GLdouble det;
  255.    GLdouble tmp[16]; /* Allow out == in. */
  256.  
  257.    if( m41 != 0. || m42 != 0. || m43 != 0. || m44 != 1. ) {
  258.       invert_matrix_general(m, out);
  259.       return;
  260.    }
  261.  
  262.    /* Inverse = adjoint / det. (See linear algebra texts.)*/
  263.  
  264.    tmp[0]= m22 * m33 - m23 * m32;
  265.    tmp[1]= m23 * m31 - m21 * m33;
  266.    tmp[2]= m21 * m32 - m22 * m31;
  267.  
  268.    /* Compute determinant as early as possible using these cofactors. */
  269.    det= m11 * tmp[0] + m12 * tmp[1] + m13 * tmp[2];
  270.  
  271.    /* Run singularity test. */
  272.    if (det == 0.0) {
  273.       /* printf("invert_matrix: Warning: Singular matrix.\n"); */
  274.       MEMCPY( out, Identity, 16*sizeof(GLdouble) );
  275.    }
  276.    else {
  277.       GLdouble d12, d13, d23, d24, d34, d41;
  278.       register GLdouble im11, im12, im13, im14;
  279.  
  280.       det= 1. / det;
  281.  
  282.       /* Compute rest of inverse. */
  283.       tmp[0] *= det;
  284.       tmp[1] *= det;
  285.       tmp[2] *= det;
  286.       tmp[3]  = 0.;
  287.  
  288.       im11= m11 * det;
  289.       im12= m12 * det;
  290.       im13= m13 * det;
  291.       im14= m14 * det;
  292.       tmp[4] = im13 * m32 - im12 * m33;
  293.       tmp[5] = im11 * m33 - im13 * m31;
  294.       tmp[6] = im12 * m31 - im11 * m32;
  295.       tmp[7] = 0.;
  296.  
  297.       /* Pre-compute 2x2 dets for first two rows when computing */
  298.       /* cofactors of last two rows. */
  299.       d12 = im11*m22 - m21*im12;
  300.       d13 = im11*m23 - m21*im13;
  301.       d23 = im12*m23 - m22*im13;
  302.       d24 = im12*m24